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Abstract 

We present a path integral method to derive closed-form solutions for option prices in a stochastic 
volatility model. The method is explained in detail for the pricing of a plain vanilla option. The 
flexibility of our approach is demonstrated by extending the realm of closed-form option price 
formulas to the case where both the volatility and interest rates are stochastic. This flexibility 
is promising for the treatment of exotic options. Our new analytical formulas are tested with 
numerical Monte Carlo simulations. 
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I. INTRODUCTION 



Since the seminal work of Black and Scholes [H, 0], who drew an analogy between the ran- 
dom motion of microscopic particles and the unpredictable evolution of stock prices, methods 
from theoretical physics have proved very useful for pricing various financial derivative prod- 
ucts jl, 0, @]. The pricing of derivative products is based on a model for the evolution of the 
probability function of the underlying asset. In order for a model to describe the economic 
reality accurately, a sufficiently general evolution for the probability distribution has to be 
allowed for. Nevertheless, the simple diffusion model of Black and Scholes (BS) is still widely 
used. Much of its success is due to the availability of closed-form analytical pricing formulas 
for many types of derivatives 

It is known for a long time that the BS model is only a crude approximation to the 
economic reality and that its assumptions are violated in actual markets. Perhaps the most 
illustrative violation is that the volatility implied from traded vanilla options, the implied 
volatility is not constant across strikes and maturities. Examples of models that tackle such 
violations are local volatility processes 0, 0], jump processes 0, Levy processes 0] and 
stochastic volatility models [lOj. A stochastic volatility model that has been particularly 
successful at explaining the implied volatility smile in e quity and foreign exchange markets 
is the Heston model In his seminal paper, Heston derived a closed form solution 
for the price of a vanilla option, which enables a quick and reliable calibration to market 
prices, especially for liquidly traded vanilla options with maturities between 2 months and 2 
years [l2| . Contrary to the Black-Scholes model, to date in the Heston model no closed-form 
analytic formulas have been found for exotic options (for recent results see (HI). Since no 
such formulas are available in the literature for any but the simplest payoffs, often costly 



numerical techniques have to be used (see [14( and references therein) 



The original mathematical solution of the option pricing problem was formulated within 
the framework of partial differential equations, but an equivalent description with path 
integral methods was developed in the pioneering work by Linetsky [l5| and Dash [l6|, 17 . 
They showed that path dependent exotic options can be straightforwardly priced with the 
path integral method. This should be intuitively clear: in the path integral formalism, a 
probability is assigned to every evolution path of the asset. In the formulation with partial 
differential equations, such quantities are typically difficult to access. 

Path integral methods have also been used in the pricing of options within stochastic 
volatility models [l8|, Ji| and in the related problem of non-Gaussian diffusion [20] (at the 
end of Sec. Ill Al we come back to this connection), but to the best of our knowledge no 



explicit option pricing formula as cheap to evaluate as Stein and Stein's [2JJ or Heston's 
formulas [ll[ have yet been derived using path integrals. We will show in the present paper 
how to carry out this task for the Heston model. The result we thereby obtain corresponds 
to the existing result for which the calibration and correspondence to market data has 
already been investigated see for example j22|, 23, 24 . 12R 26l . 27]. For a thorough discussion 
on when which approach should be used we refer to [28( and references herein. 

It is also known that there are still important features of asset price distributions which are 
absent in the Heston model for example: empirical studies of time series provide evidence of 



the long time memory of volatility |29j, ]30J • Since models containing a memory effect through 



retarded interaction, for example in the context of polarons, |3jj, have been solved within 
a path integral framework, we think our method can prove to be useful in more realistic 
models for the market also. 
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The full power of the path integral method becomes clear, when we exploit its flexibility 
by calculating the price of an option in a setting where not only the volatility but also the 
interest rate is stochastic and follows the widely used CIR model j32|, 33, 34, 3^]. To the 



best of our knowledge, no exact closed-form formula for this problem is available. Therefore, 
we have checked our formulas against a Monte Carlo simulation. 

The plan of the paper is as follows. In Sec. Ill Al we outline our model, which is the one 
introduced by Heston. Extensions of our method to different models are however straight- 
forward. Further in this section we derive a closed-form solution for the time evolution of the 
asset price. In Sec. Ill Bl we present a closed-form pricing formula for plain vanilla options 
which only involves one numerical integration of a compilation of elementary functions. In 
Sec. Ill we will extend the Heston model to include stochastic interest rate, in Sec. Ill A 
we present a closed-form solution for the vanilla option price which still contains only on 
numerical integration of a compilation of elementary functions. In Sec. IIII Bl we test this 
result with a Monte Carlo method and discuss the relevance of including stochastic interest 
rate. Conclusions are drawn in Sec. IIVI 

II. STANDARD HESTON MODEL 

A. The model and its path integral representation 

We will concentrate on assets following a diffusion process described by the following two 



equations introduced by Heston 11 



dS = fi Sdt + S\fvdwi, (1) 
dv = kq (6o — v) dt + <j\[v (pdwi + a/1 — p 2 dw7 2 j • (2) 

Here S is the asset price, /xq is a constant drift factor, v is the variance of the asset, kq is 
the spring constant of the force that attracts the variance to its mean reversion level 6o(also 
called the mean reversion speed), a is the volatility of the variance, and Wi and are 
independent Wiener processes with unit variance and zero mean. The asset price follows a 
Black-Scholes process [l[, whereas the volatility obeys a Cox- Ingersoll- Ross process (32|. 

There are two general approaches to determine the price of an option in a path integral 
context. One could, based upon equations (pQ), (j2D determine the probability distribution for 
the asset price at the strike time T conditional on the values of the asset and the variance 
at the present time Ps (St, vt \ So, vo). The expectation value of the option price at time T 
can be calculated by integrating the gain you make with a certain outcome of St multiplied 
by the probability of obtaining that outcome Ps (St, vt \ So, v o) over all possible values of 
St- To obtain the present value of the price one then discounts this expectation value with 
the risk free interest rate r. For a European call option this can be written as: 

+oo 

C = e~ rT J dSTdvT max [St — K, 0] Ps (St, vt \ So, v o) . (3) 

— oo 

We will refer to this approach as the " asset propagation approach" since Ps is the propagator 
for a distribution of asset prices (and volatilities). 
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The other approach focuses on the option price rather than the asset evolution, as will 
be referred to as the "option propagation approach" . In his paper [6], Heston discusses the 
subtle differences between the asset point of view and the option price point of view, and 
this discussion is also relevant to the present path-integral framework. Heston motivates 
that the time evolution of the option price U(S,v,t) is governed by the following partial 
differential equation (pde): 

d JL = - rS d JL + rU - {Ko [«. - „] - An} «> - >S " WS^L - ±<*,™ (4) 



dt OS J dv 2 OS 2 r dSdv 2 dv 



where A is a parameter introduced [U] on the basis of no-arbritage arguments and setting 
up a risk-free portfolio. If one makes the substitution U = e rt V one obtains the following 
pde for V as a function of the asset price and the volatility: 

dV dV dV 1 2 d 2 V d 2 V 1 d 2 V 

~dt = ~ rS dS ~ {K ° [6 ° - V] - XV} ^- 2 VS ~ pavS dSd~v ~ 2° (5) 

Based on this pde, one can find a kernel Py that propagates a given final distribution 
V (St, v, T) backwards to the present value V(5o, v , 0) of the option. Since the value of the 
option at the final time T is known, V (St, v, T) = e~ rT U (St, v, T) = e~ rT max [St — K,0], 
the value of the option now is obtained through 

+ 0O 

C = e~ rT J oIStoIvt max [St — K, 0] Py (St, vt \ S , v ) , (6) 

Furthermore the pde (j3J) is equal to the Kolmogorov backward equation corresponding to 
the following system of stochastic differential equations 

dS = rSdt + S\/vdwi, (7) 
dv = {no [9o — v] — Xv} dt + a\pv {^pdw\ + y/l — p 2 dw2^ . (8) 

This means that both approaches can be dealt with simultaneously by considering a 
generalized stochastic process: 

dS = /iSdt + Sy/vdwi, (9) 
dv = k(9 — v) dt + o\fv i^pdw\ + \/l — p 2 dw2j ■ (10) 

and calculating its transition probability P (St,vt \ So,vq). The "asset propagation" ap- 
proach ([3]) can then be retained by simply replacing /i, k and 9 by /x , ^ and 9 and the 
"option propagation" approach (jH]) by replacing //, k and 9 by r, k +X and k 9 / (k + A). 
The pricing formula for the European call is the same as ([3]) where this time the transition 
probability P (St,vt \ S ,vq) is the one corresponding to (JHJ), (TTOT) : 

P s (S T ,v T I ^o,^o) = P(S T ,v T I S ,v ) l^,^^ 
P v (S T , v t I S , v ) = P (S T , v T \ S , v ) \^ K=KO+X:9=KodQ/{KQ+x) 

We will calculate the transition density P for the general stochastic process 0, tJTU]) . 
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For later convenience we make the following substitutions: 

S 



x = In 



Sr 



~ fit, 



(12) 



x is called the logreturn and z is the volatility of the asset price. After these substitutions, 
Eq. (j2J) becomes: 

~;2 



dx = ——dt + zdw\, 



dz 

The substitution 







z 






— K — 






2_ 



dt + ^ {^pdwi + y/l — p 2 dw^ 



(13) 
(14) 



y(t) = x (t) - - (z 2 (t) - nBt) , 
a 

leads to two uncorrelated equations: 

( — K — - J Z 2 (it + Zv/l — pPdlVi, 

\a 2) 



dy 

dz 



2z 





z 


(-?) 


— K— 




2_ 



(it H — dwo- 
2 



(15) 
(16) 



We will assume that the initial volatility z(i = 0) = z$ is known [12|. The probability that 
1/ has the value and z the value zt at a later time T will be denoted as P (yx, zt \ yo, Zq). 
The advantage of transforming to these variables is that dw\ and dw2 are uncorrelated, so 
that the following expression holds for P (y?, zt \ yo,zo)' 

P(y T ,z T | y ,z ) = J VyVzexp I- J {C Q [y(t), z(t)\ + C C m{z{t)\} dt \ . (17) 



Where the quadratic Lagrangian Cq (y(t),z(t)) equals 

£ Q (^),*(t))= 2 ^_ p2) 
and the Lagrangian corresponding to the CIR process, Ccm[z(i)], is given by (36| : 









v - 











(7- 



- I k9 — ) — KZ 



1 

4^2 



K 



(18) 



(19) 



The first step in the evaluation of Eq. ( JT71) is the integration over all y-paths. Because 
the action is quadratic in y this integration can be done analytically and yields 



P (VT, Z T I yo, Zq) 



Vz(t) 



1 



exp 



v /27T^(l_p2) - | (l_ p 2) 

(z/r - i/o) 2 



a 2 



2(1-P 2 

T 

dt Ccm[z{t)\ \ . 

o 



z 2 



2z 2 (1 - p 2 ) 



(20) 
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Note that the probability to arrive in (yT, zr) only depends on the average value of the 
volatility along the path z(t): z 2 = j Q T z 2 (t)dt, in agreement with Ref. [21]. However, this 
average value appears in the denominator of the third term, and to perform the functional 
integral one needs to bring this into the numerator. This is achieved by rewriting part of 
the expression (|2T)j) as follows: 



1 



^2uz 2 (1 - p 2 ) 



exp 



(vt - yoY 

2z 2 (1 - p 2 



dk 
Ztt 



i (vt - yo) k 



Jz 2 dt(l-p 2 



-k 



(21) 

Combining Eqns. f[2"Uj) and f[2"Tj) and making the substitution k — I + i (i_ p 2) the transition 
probability becomes 



f +oc dl f 
P (Vt, Zt I yo, z Q ) = J — exp [i (y T - y ) 1} J Vzit) 



(22) 



x exp f- J dt lc CIR [z(t)} + ^z 2 



(1-p 2 ) l 2 + 2U[^-t 



The path integral over the CIR action is formally equivalent to the exactly solvable radial 
harmonic oscillator J3?J and, fortunately, adding terms proportional to z 2 to the action 
does not spoil this equivalence. The full path integral over z(t) can be carried out without 
approximations with the following result: 



P{yr,z T | yo,z Q ) = — exp 

Z7T 



K 2 Q, 



K0 1 



2f) 



K ( 2 
2 i Z T 



!) 



+oo 



X 



exp [i (y T -y )l] y/z z T 



a 2 sinh (uT) 



x exp 



2lo 



(z + z^) coth (u>T) 



AuJZqZt 

cr 2 sinh (uT) 



dl. (23) 



where 



K 



+ ilp) +1(1 -i). 



(24) 



a 

2 v \a 

is the /-dependent frequency of the radial harmonic oscillator that corresponds to the CIR 
Lagrangian ffT9]) . After transforming back to the x variable we see that also the integral over 
the final value zt can be done analytically (see e.g. j38|), yielding the marginal probability 
distribution V (xt \ 0, z ) = f_ °° dzrP (xt, zt \ 0, z ) (written in the original variable v) as 
a simple Fourier integral: 

1 r k 

V (x T I 0, v Q ) = — exp — (kQT + v ) 



2tt 

-oo 



X 



t k0 



exp < % 



x T + - (y + k6T) 
a 



2oo 



a 2 sinh (ojT) 



[cosh (uT) - N] v \ dl 



(25) 
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where N is: 



cosh (uT) + 2b ( K + ilpo) sinh (ujT) ^ ^ 

Note the similarity of the expression (125]) with the result obtained in Ref. 2cj, derived for 
a general stochastic process with non-Gaussian noise. 



B. Pricing of plain vanilla options 

From now on we follow the option propagation approach and set p equal to r. The price 
of a call option C with expiration date T and strike K when the transition probability is 
known is given by Eq. Q. Writing this formula in the x variable and thereby inserting the 
result (|25p for the transition probability results in: 



C 



-rT 



+00 

J dxx max [S exp — K, 0] V (xt I 0, v ) 



(27) 



where the risk free interest rate was restored and denoted by r. Now there are still two 
numerical integrations that have to be done. Following the derivation outlined in Ref. [2(| 
we can rewrite expression (1271) so that only one numerical integration remains: 



C 



So - e~ rT K 



+ i 



1 ' GXP 



/ 1 —a + x e — rT) I H — -a 

a / o l . 



(28) 



with 



x 



So exp ( a ) — e r K exp (T) 



-S + e-M|, 



x 



V 

M 

e 

T 



In 



K 



So, 

a = v + k9T, 



cosh (z/T) + — (k + ilpa — pa) sinh (uT) 



Ivvo 



a 2 sinh (z/T) 

2ujv 
a 2 sinh (ujT) 



[M - cosh (z/T)] + —k6 log M, 



[N - cosh {uT)\ + ^kB log N. 



(29a) 
(29b) 

(29c) 
(29d) 
(29e) 
(29f) 



and u) defined as before (jz4[) . We have tested this result against the formula stated in Ref. 
11] . This confirmed the correctness of formula ( l28l) . Now we are confident to explore new 
grounds with our method in the following section. 
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III. STOCHASTIC INTEREST RATE 



A. Derivation of the option price 

In the previous section we assumed the interest rate to be constant. Here we allow the 
interest rate to change in time, r(£). Applying Black and Scholes' no-arbitrage argument on 
Heston's risk-free portfolio motivation for the evolution of the option price, we again obtain 
the partial differential equation (jSJ) with r(t) rather than a constant r : 

dV dV dV 1 2 d 2 V d 2 V 1 d 2 V 

■at = - r{t)s os - {Ko [9 ° - v] ~ Xv} ^- 2 vS ds^ ~ pavS dSd- v ~ 2° v w (30) 

For a given function r(i) this leads to a kernel Py [St, vt \ So, v \ r(t)] so that the option 
price becomes 

+oo 

C[r(t)} = J dS T dv T max [S T - K, 0] e^ r( ' )d 'P y [S T ,v T \ S ,v \ r(t)} . (31) 

— oo 

Note that the option price is now a functional of the given time evolution of the interest 
rate r(t). As in the previous section, it is convenient to introduce new integration variables 

V (t) = In (jQ - P - [z 2 (t) - K9t] , (32) 

z(t) = yfitf). (33) 

In the path- integral treatment, the kernel can be written as a sum over all possible realiza- 
tions of y(t) and z(t), weighed by the action functional of the system: 

+oo 

C[r{t)] = J dx T dv T max [e XT - K, 0] e~fo r[t)dt 

— oo 

x J VyVz exp ^- J {C Q [y(t), z(t), r(t)} + C C m[z(t)}} dt j , (34) 
where Cq is the quadratic Lagrangian (|T8|) 



C Q (y(t),z(t)) ' 



2z 2 (1 - p 2 



m-r(t)-[^-^z 2 (t)T , (35) 



and Ccir is the CIR Lagrangian. Of course, we cannot know what particular realization of 
the interest rate r(t) will appear in the future. We assume the interest rate to follow a CIR 
process which is uncorrelated from the other two stochastic processes, 

dr = ) dt + a r y/rdw 3 . (36) 

The value for the option price then needs to be averaged over the realization of r(t) in this 
CIR process. Where the calculation of the expectation value of such a functional might 
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become cumbersome with conventional probabilistic techniques, it can be evaluated very 
elegantly with the Feynman-Kac formula: 



- J C OIR [r(t)}dt 



(37) 



where Cqir is the Lagrangian for the CIR process. The final result can be expressed with 
a modified propagator P(St,vt,tt \ So,vo,r Q ) as 



-l-oo 

C = J dSxdvTdr T max [St — K, 0] P(St, vt, r T \ S , v , r ) 



(38) 



with 



P(S T ,VT,r T \So, v ,r ) = / VyVzVr e 



-So'r{t)dt 



x exp 



t T \ 

-J {C Q [y(t),z(t),r(t)}+CciR[z(t)}+Ccm[r(t)}}dt\ . (39) 



The stochastic interest rate makes the vanilla price dependent on the specific path followed 
by the interest rate. This part of the payoff has been taken into the calculation of the 
propagator, where it is analytically tractable, and no longer appears explicitly in the ex- 
pression ([38"]) for the option price. Herein lies the strength of the path-integral approach, 
to price path-dependent options. With a stochastic interest rate the European vanilla op- 
tion becomes dependent on the entire path of the interest rate and is still solved in a very 
straightforward way. This is promising for more general option types, such as the barrier 
and Asian options that we are currently investigating. 

A useful substitution to perform the functional integrations is 



Mt) = y(t)- I r(t')dt'. 



(40) 



As was the case for the Lagrangian corresponding to the volatility, the Lagrangian cor- 
responding to the interest rate process will also be formally equivalent to the Lagrangian 
corresponding to a radial harmonic oscillator; furthermore the addition of another term 
quadratic in $i stemming from the discount factor doesn't spoil the correspondence. The 
result reads as follows: 



+ i 



S - if exp ( -^a r + T r (0) 



err. 



K exp 



T r (0) + —a r 



- So + exp 



-a + x t 



S exp (-^a + 6 + e r ) - if exp (T + T r ) } 



a 

dl 
2k' 



I + —a 

a 2 



(41) 
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To make it surveyable, we introduced the following notations 



a r = tq + K r 9 r T, 



io r (I) 



:) = ^| + 2( ii + l), 



-i -1 



M r = cosh (v r T) + 



2u r 



sinh (v r T) 



[M r - cosh (v r T)) + 2^f log M, 



T r (I) 



of sinh (u r T) 
2u r (I) r 



{N r (I) - cosh [u r (I) T]} + 2-^ \ogN r (I) . 



r 



of sinh [u r (I) T] 



r 



These notations reflect the extension to the case of stochastic interest rate (symbols with 
subscript r) of the corresponding quantities in the Heston model (equations (I29ap - fl29fj) ). 
Notice the resemblance with formula (1281) . Formula (HT1) still contains just one numerical 
integration with an integrand composed out of elementary functions. To the best of our 
knowledge, only approximate analytical formulae are available when both the volatility 
and interest rate are stochastic [39(]. Because of the lack of alternative exact analytical 
expressions, we have checked the validity of our formula (|41|) against numerical Monte Carlo 
simulations. Our Monte Carlo method is outlined below. 

First notice that substitutions (140!) transform the x- variable into a variable x, independent 
of the interest rate by subtracting the time averaged interest rate f: x = x — f. This results 
in the same equation as in the constant interest rate situation, Eq. (|T3l) . Also the discount 
factor only contains f. This means that the knowledge of the probability distribution f is 
sufficient to calculate the price by means of the formula (I28p derived in the constant interest 
rate setting. So the Monte Carlo scheme used is the following: first values for f are simulated 
and used to calculate the option price for these values, next the price is averaged over all the 
simulations. A value for f is simulated as follows: time is discretized in little time steps A, we 
sample a path for r and integrate along this path. To calculate the probability distribution 
for f, we used the result that the stochastic time increment of a CIR variable over a small 
time step At follows a non-central y 2 distribution [3]. The probability distribution of the 
average interest rate f is then simulated by generating many r-paths in discretized time. As 
shown in Fig. [H the agreement between the analytical (thick full line) and numerical option 
prices is excellent. 

In this section the option propagation approach was followed from the beginning. In this 
setting it is necessary to make a choice between the two approaches from the start because 
in the asset propagation approach one would actually have to introduce a stochastic process 
for the drift /io instead of for the interest rate. That these two should follow the same 
stochastic process is not clear. Since the option propagation approach is the most common 
one anyway we followed this approach. If one does want to introduce a stochastic process 
for the drift //q this poses no problem and the derivation of an option price in this setting 
would be completely similar. 
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0.6 0.8 1 1.2 1.4 0.6 0.8 1 1.2 1.4 0.6 0.8 1 1.2 1.4 



K/S Q K/S K/S 

FIG. 1: (color online, two column wide) This figure shows the result of different pricing formulas 
from which the Black-Scholes result (with interest rate r = 9 r ) has been subtracted. Since we are 
not considering a specific asset, the option price could be stated in any currency, therefore the 
deviation is given in arbitrary units. The thick (red) curve shows our analytical results for the 
model with both stochastic interest rate r and stochastic volatility. The crosses represent results 
from a Monte Carlo simulation of our model, confirming the analytical formula. The blue dotted 
curve and the blue dash-dotted curve show the results for the Heston model with constant interest 
rate r = r(0) and r = 9 r , respectively. The dashed curve shows the results for a Black-Scholes 
model with r = r(0). The following parameter values were used for all three panels: k = 1, a = 0.2, 
9 = 0.04, v = 0.04, T = 1, S = 100, K r = 1.8; a r = 0.1; 6 r = 0.03; r = 0.035. The correlation 
coefficient is for panel (a) p = —0.5, (b) p = and (c) p = 0.5. 



B. Results and discussion. 

In the current treatment, we have two layers of generalization as compared to the Black- 
Scholes result. First, the volatility appearing in the Black-Scholes process is stochastic - 
this leads to the Heston model. Second, the interest rate of the Black-Scholes model is 
also stochastic - leading to our present results. In this paragraph, we argue that both 
improvements can have an equally important effect on the option price. 

This is illustrated in Fig. [TJ where the different approaches are compared. Let's start 
with the most complete model, where both interest rate and volatility are stochastic. The 
resulting option price, Eq. f|4Tl) . for this model is shown as a thick red curve. The result from 
the closed-form expression agrees well with the Monte Carlo simulation, shown as crosses. 

Now we strip off one layer of complexity, and fix the interest rate r - it is no longer a 
stochastic variable. Then we obtain the Heston model as an 'approximation' to a stochastic 
interest rate world. The question poses itself of which fixed interest rate to use, if we want 
to make the comparison. Two choices are shown in Fig. [TJ r = r(0) and r = 6 r . The 
former choice (dotted blue curves) sets the Heston interest rate equal to the interest rate 
at time 0, whereas the latter choice (dash-dotted curves) sets the Heston interest rate equal 
to the mean reversion level 6 r . For the parameter values used in Fig. [TJ the most complete 
result lies between the two Heston 'approximations', but this is not necessarily so. Fig. [2] 
shows that for some choices of other (realistic) parameters, the full result can lie outside 
both Heston approximations. Nevertheless, as k becomes very large, the stochastic interest 
rate will be drawn very tightly to the mean reversion rate 8 r , and one expects the full result 
to be near the Heston approximation with r = 6 r . When k is very small, the stochastic 
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0.5 



Heston, r = 9, 




-0.2 - 

0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 
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FIG. 2: (color online). As in Fig. 1, the result of different pricing formulas from which the Black- 
Scholes result (with interest rate r = 9 r , thin black line) has been subtracted, is shown. The 
following parameter values were used: n = 1, a = 0.2, 9 = 0.04, vq = 0.04, T = 1, So = 100, 
K r = 0.5, a r = 0.3, 9 T = 0.03, r = 0.035, p = 0. 



interest rate will not be drawn quickly towards 9 r so that when also a r is small, the full 
results will be near the Heston approximation with r = r(0). 

Next, we strip off the second layer of approximation, and also fix the volatility. This 
results in the familiar Black-Scholes model as the crudest approximation to our system. 
Now a second choice has to be made: which value of the volatility to use. Here, we take the 
stochastic volatility at time zero to be equal to the mean reversion level of the volatility CIR 
process, so that the ambiguity of choice is avoided. The choice for what interest rate to use, 
however, remains. In Fig. [T], we show the Black-Scholes results with r = r(0) (dashed line) 
and r = 6 r (full line). We have plotted all the results relative to the Black-Scholes result 
with r = 9 r to emphasize the differences rather than the absolute magnitude of the prices 
(for this reason, the r = 8 r Black-Scholes result is the baseline of the plots). The difference 
between the three panels of Fig. [T] is the value of the correlation between asset price and 
volatility. 

From Figs. [T] and [2J it is clear that both levels of approximation (keeping the volatility 
constant and keeping the interest rate constant) have an equally large effect on the option 
price. Even within the Heston framework, the choice of what value to use for the interest 
rate is seen to influence the price considerably. Choosing a different interest rate, or keeping 
the interest rate as a stochastic variable, leads to a price correction that is as large as the 
price correction obtained by going from the Black-Scholes to the Heston model. This result 
emphasizes the importance of a correct treatment of the interest rate in pricing models (This 
also depends strongly on the length of the lifetime of the option). 

Finally we must remark that the price differences when working within the standard 
Heston model or within the extended one can be influenced by the calibration method. 
For Figures. [1] and [2] we used the same parameters for the volatility process both in the 
standard model and in the extended one, parameter values for the interest rate process are 
calibrated separately. Literature shows that the parameter values for the volatility process 
(see for example j24| and [ll[) and the interest rate process (see for example [35| and (33j ) 
can attain values in a broad range containing the values we chose to produce Fig. [1] and Fig. 
[2j However if the parameter values obtained for the interest rate process are used in formula 
(jUJ) to calibrate the remaining parameter values for the volatility process one might get 
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different results. We can not exclude that this calibration approach would lead to smaller 
price differences between the two approaches. However such a calibration is a research area 
on its own and is outside the scope of this article. 

IV. CONCLUSIONS 

We have developed a path-integral method to derive closed-form analytical formulas for 
the asset price distribution in the Heston stochastic volatility model. Closed-form formulas 
are obtained for the logreturn of the derivative and the vanilla option price. The presented 
results correspond to the known semi-analytic results obtained from solving the partial 
differential equation [ll[ by standard techniques. 

The flexibility of our approach is demonstrated by extending the results to the case where 
the interest rate is a stochastic variable as well, and follows a CIR process. For this case, 
to the best of our knowledge, no exact analytical solutions have been derived before. We 
have checked our semi-analytical results for the model with both stochastic volatility and 
stochastic interest rate against a Monte-Carlo simulation. The quantitative analysis shows 
that the effect of stochastic interest rate on the Heston model can be as large as the effect 
of the stochastic volatility on the Black-Scholes model. However we did not perform a 
full calibration, which might influence the results. Finally, the analogy between stochastic 
interest rate models and path dependent options makes our method promising for the pricing 
of exotic derivative products. 
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